Exploring Valleys of Aging Systems: The Spin Glass Case 
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We present a statistical method for complex energy landscape exploration which provides infor- 
mation on the metastable states — or valleys — actually explored by an unperturbed aging process 
following a quench. Energy fluctuations of record size are identified as the events which move the 
system from one valley to the next. This allows for a semi-analytical description in terms of log- 
Poisson statistics, whose main features are briefly explained. The bulk of the paper is devoted 
to thorough investigations of Ising spin glasses with Gaussian interactions of both short and long 
range, a well established paradigm for glassy dynamics. Simple scaling expressions with universal 
exponents for (a) barrier energies, (b) energy minima, and (c) the Hamming distance as a function 
of the valley index are found. The distribution of residence time inside valleys entered at age is 
investigated, along with the distribution of time at which the global minimum inside a valley is hit. 
Finally, the correlations between the minima of the landscape are presented. 

The results fit well into the framework of available knowledge about spin glass aging. At the same 
time they support a novel interpretation of thermal relaxation in complex landscapes with multiple 
metastable states. The marginal stability of the attractors selected is emphasized and explained in 
terms of geometrical properties of the landscape. 

PACS numbers: 05.40.-a, 65.60.+a, 75.10.Nr 



I. INTRODUCTION 

Physical properties of glassy systems quenched from a 
high temperature slowly change with the time or age 
elapsed since the quench. For any well below the true 
equilibration time, the dynamics deceptively appears to 
be stationary when observed on time scales shorter than 
tiu. Experimental and numerical evidence for the pres- 
ence of a quasi-stationary fluctuations regime for t <^tw 
followed by non-equilibrium drift for t ^ t^, stems e.g. 
from measurements of conjugate linear response and au- 
tocorrelation functions [E S 0] which obey, respectively 
violate, the fluctuation-dissipation theorem in the two 
regimes. In general, the (apparent) age of the system 
can be deduced from the decrease in the rate of change 
of macroscopic averages. This apparent age can however 
be 'reset' to an earlier value by applying a perturbation 
of short duration, as e.g. a temperature pulse, which thus 
rejuvenates the system. 

Non-equilibrium memory effects such as aging and re- 
juvenation were first noticed and studied in a spin glass 
context 0, 0, 0,13, but are now observed in a variety of 
glassy systems ITol [ill IT^ . Expanding a previous 

brief exposition |l3j| . we describe and test a statistical 
approach to landscape explorations designed to find the 
generic properties of the energy landscape which produce 
these effects. 

Our main point is that the dynamical events mark- 
ing the transition between the quasi-equilibrium and the 
off-equilibrium dynamical regimes are triggered by the 
attainment of energy values of record magnitude. This 
immediately allows a description of the non-equilibrium 
dynamics in terms of a log-Poisson process 14, 15], i.e. 
a stochastic process which is homogeneous when viewed 
on a logarithmic time scale. 



The sequel is organized as follows; the next section in- 
troduces the landscape exploration method. Section UTTI 
briefly explains the relevant properties of the log-Poisson 
statistics used throughout the paper as a semi-analytical 
description of non-equilibrium dynamics. Section 
presents results of an extensive application of the method 
to spin glass systems, demonstrating the viability of the 
method and the usefulness of the log-Poisson description. 
Section lVll Duts the results in a broader perspective, with 
special reference to coarse-grained mesoscopic models of 
conflguration space. In particular, we discuss the con- 
nection between aging in thermalizing systems and in 
dissipative driven system and biological evolution. Sec- 
tion is a summary and an outlook. 

As the nature of true equilibrium and the final stages 
of thermal relaxation are only weakly related to the pre- 
sented data, no theoretical consideration is given to these 
aspects in the paper. 



II. THE METHOD 

In models of driven dissipative systems with multiple 
attractors 0, 0, , marginally metastable attractors 
with an a priori negligible statistical weight are never- 
theless those typically selected by the dynamics. In such 
systems, this mechanism underlies memory and rejuve- 
nation effects analogous to those observed in the ther- 
malization of e.g. spin glasses after a quench. One 
can speculate that similar mechanisms could generally 
be present in glassy systems with an extensive number 
of metastable attractors. However, the issue of attractor 
selection is not explicitly considered in widely used land- 
scape exploration methods such as the Stillinger- Weber 
approach El El S IH HI , which study a set of local 
energy minima (inherent states) generated by quenches. 
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The same applies to exhaustive landscape exploration 
techniques [2^ 0; IHl and studies of the real space mor- 
phology of low-energy excitations by techniques requiring 
quenches |2^, genetic algorithms '27"| and energy mini- 
mization of excitations of fixed volume ,28i . .29,] . 

The ability of small external perturbations to induce 
strong rejuvenation and memory effects in complex dy- 
namics strongly suggests that any probe introducing ex- 
traneous elements in the dynamical evolution might at 
the same time yield a biased picture of the energy land- 
scape. In other words: are the attractors identified also 
those which would be selected by e.g. the unperturbed 
thermalization process after a deep quench or any other 
dynamical evolution of interest? Analyzing the regions of 
state space surrounding intrinsic states provides valuable 
information about the quasi-equilibrium (fluctuation) dy- 
namics in the energy landscape, but it docs not tell the 
whole story. The question of how to properly describe 
the non-equilibrium process of 'selecting' the metastable 
states is still open. This question motivates the present 
approach which is based solely on statistical information 
collected during an undisturbed aging process. 

Conventionally, a valley is a connected neighborhood 
of configuration space which supports a state of approxi- 
mate thermal equilibrium centered on a local energy min- 
imum. During the time a trajectory 'resides' in a valley, 
the energy and other physical quantities fluctuate around 
a fixed average, and the state of lowest energy is often 
revisited; the dynamics has a recurrent character. Non- 
equilibrium events, henceforth 'quakes', move the system 
irreversibly from the neighborhood of the initial local en- 
ergy minima of high value and into progressively deeper 
valleys. A sequence of such events seldom or never revis- 
its the same configurations and has a transient character. 

As widely recognized, the lack of time translational 
invariance in aging systems stems from the dynamical 
in-equivalence of the valleys visited. Consider therefore 
a landscape with multiple valleys of varying degrees of 
metastability, or depth. On time scales larger than the 
residence time of the deepest (i.e. most stable) valley seen 
up to the age t, all valleys shallower than this valley are, 
by definition, unstable and hence irrelevant for the non- 
equilibrium dynamics. The interesting dynamical objects 
arc thus the valleys deeper than the deepest valley seen 
up to time t. This points to energy records as prospective 
markers of the non-equilibrium events. 

We shall use the term 'energy barrier' to denote the 
difference between the energy of the current state and the 
lowest seen or best-so-far energy minimum. The lowest 
energy minimum value and the highest energy barrier 
observed up to time t will identify the valleys as they 
successively appear in the landscape. 

Our classification procedure of the undisturbed dy- 
namics is as follows: We save the minima and barriers of 
record values encountered and the times at which they 
occur. We stipulate that a valley is entered at time t if 
the barrier record achieved at time t happens to be the 
last one prior to the attainment of a energy minimum 



record. To leave a valley, a barrier record must again 
be followed by a record in the lowest energy. Whenever 
several minima records are found between two barrier 
records, we only keep the latest, and therefore deepest, 
minimum. 

In short, we operationally identify the valleys encoun- 
tered by a series of minima records with at least one bar- 
rier event between them. Note that this selection proce- 
dure must be performed retrospectively, since it is impos- 
sible to know 'on the fly' whether a new valley has been 
entered or not. We also stress that the barrier records 
will not necessarily coincide with the lowest barrier sep- 
arating two consecutive energy minima records, as the 
dynamics, due to entropic effects, is not likely to follow 
the path of lowest energy, an effect noticed by Wevers et 
al. | 30l | in the landscape of metastable ionic compounds. 

The discovery of a new record in low energy is a non- 
equilibrium event. However, by no means does it im- 
ply that the internal dynamics in a valley is entirely 
equilibrium-like. Several sub-valleys are typically ex- 
plored before the energy minimum is encountered which 
eventually remains as the lowest state within the val- 
ley. Only then does the dynamics acquire the recurrent, 
fluctuation-like nature which is characteristic of quasi- 
equilibrium. 

Resetting the highest barrier to zero at an arbitrary 
point in time produces numerous barrier records but no 
new valleys before a lower energy value is again recorded. 
However, resetting both the record energy and barrier 
values to zero may result in a series of new records being 
registered, which describe sub-valleys within the valley 
originally explored. By repeating the simulation with 
the exact same random numbers, this procedure allows 
one to take a closer look at the internal structure of a 
valley if the resetting is done at the time of entry. 

The method presented is generally applicable, easily 
implemented, and does not add much to the total run- 
time of a simulation. If the energy landscape explored is 
simple, e.g. if it contains one large, structureless valley or 
if it is perfectly periodic, our scheme only detects a single 
valley, since degenerate and hence dynamically equiva- 
lent minima are appropriately lumped together. In other 
words, our method produces simple results when used on 
simple systems. In the following sections we show that 
non-trivial results are indeed the outcome when complex 
landscapes are explored. 



III. LOG-POISSON STATISTICS 

In this section we introduce and motivate the log- 
Poisson statistics used through the paper as an idealized 
analytical description of non-equilibrium events — such as 
the quakes which lead to new valleys. We mainly focus 
on the consequences and predictions for various quanti- 
ties of physical interest, leaving the empirical justification 
of the formalism to the next section, where the data are 
also presented. 
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A log-Poisson description of complex dynamics was 
first introduced in connection with a model of Charge 
Density Waves and later used to explain macro- 

evolutionary patterns from the fossil record 'sT, "s^]. It 
describes j33| the coarse grained dynamics of a popula- 
tion evolving in a NK landscape and there are in- 
dications that it might also apply to far more realistic 
models of evolution 35]. 

The familiar Poisson process with the time argument 
replaced by its logarithm is in short denoted log-Poisson: 



P.W = ^^-", t>i. 



(1) 



As shown in 0, [s^l , the probability that k records occur 
in a sequence of t random numbers is given by Eq. 
independently of the underlying process generating the 
numbers. Log-Poisson statistics implies that the tempo 
at which the events occur falls off as 1/t. Switching to 
logt as the independent variable gives a constant (log- 
arithmic) rate of events, and restores time homogene- 
ity. Several other interesting mathematical properties of 
log-Poisson processes are listed below together with their 
physical implications for the spin glass problem at hand. 

The probability that k 'events' occur between tw and 
tui -Ht is ITU 



a In 



ill! 



(2) 



Consider a function c(m, m + k) describing the effect of 
k events, for example the overlap between the configura- 
tions of lowest energy in valleys m and m + k. To connect 
with the thermal correlation C {t^ ,tw+t), we can utilize a 
cartoon rendering of the coarse-grained non-equilibrium 
dynamics as a one dimensional walk on a set of states cor- 
responding to the minima of the valleys, which is similar 
to e.g. one of the basic assumptions of Ref. [s^. We as- 
sume, in our terminology, that the dynamical trajectories 
of an aging system dwell at the bottom of the 'current' 
valley for a certain residence time, until a quake instan- 
taneously moves them to the bottom of the next valley. 
This coarse grained picture neglects the internal struc- 
ture of the valleys, which is questionable since the time 
spent searching for the bottom state within a valley is of 
the same order of magnitude as the full residence time, 
see e.g. Fig. [7| Nevertheless, the cartoon has the virtue 
of simplicity, and allows us to write 

C(t^,t^, + t) Prn{tw)Pk{tw,tw + t)c{m,m + k). (3) 



which is a function of t/tw if and only if c(m, m + k) is 
independent of m. This pure or full aging form has very 
recently been shown to describe aging in real spin glasses, 
if cooling rate effects through the critical temperature 
are accounted for |37j . The assumption that c(m,m -|- 
k) only depends on a single argument has been checked 



separately, and an additional dependence on the valley 
index m has been found, which entails a deviation from 
t/tu, scaling jsl] for the numerical models. 

For later convenience we finally note that, if tk marks 
the time of the k'th event in a log-Poisson process oc- 
curring after t, the 'log-waiting time' log(ife) — log(t) 
is exponentially distributed, in perfect analogy to the 
usual Poisson process. The same exponential distribu- 
tion also describes the series of independent variables 
log(ife) — log(tfe_i). This property of the log-waiting time 
distribution is very well fulfilled by our data. 



IV. MODELS AND DYNAMICS 

As glassy dynamics is insensitive to many details of 
the interactions, computational convenience is a prime 
criterion for choosing a test model. Spin glass models 
are comparatively easy to simulate and have been in- 
vestigated experimentally and numerically for more than 
twenty years. The lack of a comprehensive and coherent 
picture of their dynamics and statics furthermore endows 
them with considerable intrinsic interest. 

This section demonstrates that a simple and consistent 
geometrical picture of the spin glass energy landscape 
is obtained with our 'non-invasive' exploration method. 
Current spin glass issues are mentioned as needed, while 
a more complete discussion can be found in Section I VII 

We consider N Ising spins, where the energy of a spin 
configuration s = {si, . . . ,sn} is given by 



ijSiSj, 



(4) 



where the couplings Jy are symmetric, independent 
Gaussian variables of unit variance. To test the impor- 
tance of the topology of the system, we have simulated 
lattices with periodic boundary conditions in both 2d, Zd 
and Ad {N — i"*), where Jij are non-zero only if i and 
i are lattice neighbors. Additionally, we have placed the 
spins in a fc-regular random graph, i.e. each spin interacts 
with exactly k spins chosen at random. 

We use single spin flip dynamics coupled with a rejec- 
tionless algorithm, the Waiting Time Method Is^. The 
latter generates a sequence of moves equal in probabil- 
ity to the sequence of accepted moves in the Metropolis 
algorithm. Thus, the results of this paper can also be ob- 
tained with the standard Metropolis algorithm, albeit at 
the price of considerably longer run times. The 'intrinsic' 
and size independent time variable t used throughout cor- 
responds to the number of Monte Carlo (lattice) sweeps 
in the Metropolis as well as to the physical time of a real 
experiment. We refer to for a detailed account of the 
Waiting Time Method. In all our runs, we conventionally 
skip the data within the first 10 time units in order to let 
the system settle down from the random T = oo initial 
configuration. 
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FIG. 1: The average number of valleys nv visited as a func- 
tion of time in 3d Gaussian spin glasses with A'' = 16^ spins. 
For t > 100 a curvature is barely discernible. The insert shows 
nv for system sizes L = 6, 8, 10 and 12 (full lines), as well as 
1/ = 16 (dashed line as in the main panel). 



V. RESULTS 

Figure ^ shows the average number of valleys ny (t) 
(which is equal to the number of quakes leading from a 
valley to another) observed in 3d at various low temper- 
atures. Neglecting the earliest part of the simulations 
where the memory of the initial random spin configura- 
tion is still present, riv{t) has a logarithmic shape. As 
the ratio of the variance (t ) to riy (i) is constant and 
close to one (see Fig. 2 of ref. T5]), the statistics of val- 
leys is essentially a log-Poisson process, as expected in a 
record induced dynamics ■ As this holds true for bar- 
rier records as well (not shown) , the ratio of the number 
of barrier records to minima records is therefore constant 
on average. This non-trivial geometrical feature of the 
spin glass landscape has not previously been noted, and 
it provides a link to an analytical description of the non- 
equilibrium dynamics as a log-Poisson process 0, 0| . 

The barely perceptible curvature seen in Fig.^for large 
t decreases systematically as the system size increases at 
a fixed temperature, as shown in the insert. The sub- 
logarithmic form of nv{t) means that the likelihood that 
a record in low energy follows a barrier record decreases 
very slowly, but systematically as the system ages. Tak- 
ing into consideration that the curvature seems to vanish 
in the limit of a very large system (see Fig. 2 of ref. [T3l| . 
which shows that the slope of the straight part of riy (t) 
seems to be proportional to log(A^)), we surmise that its 
presence reflects the increasing difficulty in finding new 
low energies as the ground state is approached. 

In summary, the series of quakes moving the system 
from one valley to the next can be meaningfully ideal- 
ized into a log-Poisson process discussed in the previ- 
ous section, if one disregards the curvature of the data. 
How record events are distributed in time is insensitive 
to the properties of the stochastic process from which the 




valley index 



FIG. 2: Average energy barrier records B separating contigu- 
ous valleys in random graph spin glass landscapes are scaled 
with system size A'^ and temperature T as indicated in the 
ordinate label and plotted versus the valley index. The data 
points are based on A' = 4096 and T = 0.3, ...,0.7 as well 
as T = 0.5 and N = 1000, 16000. The same scaling form 
B ~ jv0-38yi.8 applies equally well to our 2d, 3d, and 4d data. 



records are drawn In our case Figs. Q and have a 
weak temperature dependence, as opposed to the strong 
T dependence of the underlying fluctuations. Finally, we 
note that identical results with respect to ny (i) are also 
found in 2d, id and random graphs. 

The rest of this section deals with two scaling plots 
similar to those presented in ref. [3 (Figs. |2 and as 
well as new results concerning the scaling of Hamming 
distances (Fig.|3Jl between the local minima of contiguous 
valleys. Finally, the statistics of the residence time and 
the fraction thereof spent before hitting the lowest energy 
state in the valley is analyzed, along with the correlations 
between the minima of the valleys. 



A. Scaling of barriers 

The data shown in this section are averages over many 
thousands of realizations of the couplings Jy and cover 
a wide range of low temperatures and system sizes. The 
raw data for each T and N have very little scatter, and 
the main sources of error are systematic. For example, 
our flnite runtimes of i = 10^, combined with the very 
broad distribution of residence times in the valleys, bias 
the average energy of valleys discovered late in the pro- 
cess, since only the 'faster' trajectories are able to explore 
these. The data are parametrized by the valley index i. 

In our plots the scaling of the ordinate is given in the 
corresponding labels. The value of the valley index i is 
shifted from one data set to another by up to one unit. 
This corresponds to a multiplicative shift of the age, and 
compensates for the arbitrariness of skipping data within 
the flrst 10 time units after the quench, irrespective of 
temperature and system size. 

Figure|21is a scaling plot of the energy barrier Bi which 
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on average must be surmounted in order to leave the i'th 
valley. The system is a /c-regular random graph with 
fc = 6, i.e. the same number of links per spin as in 3d. 
Bi is seen to scale in a simple way with the tempera- 
ture and the size of the system: Bi{N,T) = B,{N°-T''), 
where a — 0.38(2) and b — 1.8(1). Interestingly, the 
very same values of a and b are found when scaling the 
barrier heights in 2d, 3d and Ad lattices as well (see e.g. 
Fig. 3 in ref. 13]). So, in addition to the non trivial 
fact that a straightforward two parameter fit works very 
well, the values of the scaling exponents seem universal. 
As the next section will reveal, this universality is also 
found when looking at energy minima and the Hamming 
distance between them. 

The close to linear shape of the scaling plot shows that 
the barrier heights between valleys gradually increase. 
Again, we stress that the barriers discussed are those 
found along an unperturbed trajectory, end hence un- 
likely to be the lowest barriers between energy minima 
in contiguous valleys. 

Since the exponent a is positive, the quakes do not re- 
main localized to a finite set of spins in the macroscopic 
limit. Some information on the size and shape of the 
quakes can be inferred from the scaling law for the en- 
ergy barrier, if we assume that quakes correspond to the 
motion through the system of a generalized domain wall. 
Letting m be the number of spins typically involved in 
a quake, we write B oc for some exponent x > 0. 
Were the barrier energy the outcome of a fluctuation, i.e. 
the sum of m contributions of random sign, the behavior 
would be diffusion-like, with x = 1/2. However, since the 
barrier crossing process favors low-energy states, we ex- 
pect a lower barrier energy, i.e. x < 1/2. As m oc jSf^-^^/^ 
and m < N, we conclude that x G [0.38,0.5]. In sum- 
mary, quakes have a fractal shape with exponent 0.38/a;, 
and are more space filling than a conventional domain 
wall in d < 4. 

An explanation of the strongly non-Arrhenius T de- 
pendence of the exponent b involves entropic effects and 
the connection between the quasi-equilibrium dynamics 
in configuration space and in real space. The qualitative 
argument given below predicts that the barriers grow lin- 
early with the valley index, and that 6 = 2 which is close 
to the actual value of 6 = 1.8. The discrepancy might 
arise because we neglect that domain sizes are distributed 
quantities. 

In real space, the relevant objects are connected do- 
mains |2^ 23 of thermally equilibrated spins which 
slowly grow in a sea of frozen spins reaching, on average, 
a linear length scale £{t) on a time scale t. In the progres- 
sively longer quiescent periods between the quakes, these 
domains do not interact, and their quasi-equilibrium 
properties are therefore determined by the local density 
of states 'D{E) pertaining to the configurations accessi- 
ble to the spins constituting each domain. T^(E) has 
been investigated by means of the lid- algorithm 123] fo r 
a number of different glassy systems |24L l4ll l43l l44| , 
and has consistently been found to have a close to expo- 
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FIG. 3: The average Hamming distance H between the low- 
est energy configurations of two contiguous valleys in 3d spin 
glasses is plotted as a function of the valley index. Combi- 
nations of T = 0.5 and L = 8, 10, 13, 16, 20, 30, 40 as well as 
L — 16 and T — 0.4, ...,0.8 are included. Again, the scaling 
exponents of A'^ and T are independent of the topology of the 
networks. 

nential shape: T>{E) « exp(i?/T'o(^)). The energy scale 
To(£) is an upper bound for the temperature at which 
metastability can hold and decreases monotonically with 
the linear size of the system considered While the 
results quoted pertain to systems of fixed size and shape, 
preliminary investigations confirm that the same behav- 
ior applies to spin domains of size £(t) which grow within 
a large system. In this case To{£) becomes a slowly de- 
creasing function of time through the time dependence of 
i{t), and asymptotically approaches T from above. We 
anticipate that the marginal stability of the valley visited 
implies that a domain is typically close to its maximal 
size, i.e. To(^(i)) w T. 

Returning to Bi, we recall that energy barriers delim- 
iting valleys are extremal values in a series of 0{t) in- 
dependent outcomes in a system of age t, and that each 
attempt sees energy E with probability 

P,^{E)^exp(^-E^^y (5) 

It follows ^3 that the typical energy barrier scales as 

Scx-^^logL (6) 

Since the valley index grows linearly with logt, the linear 
dependence on the valley index seen in Fig.|21is recovered. 
Finally, since Tq « T -t- e for some small e, one obtains 
& w 2 as anticipated. 



B. Scaling of Hamming distances and energies 

The fact that contiguous valleys contain rather differ- 
ent low energy configurations is shown in Fig. |31 which 
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depicts the scaling of the Hamming distance H (num- 
ber of spins which differ in their orientation) separat- 
ing the configurations of lowest energy at the bottom of 
the two valleys in id lattices. Again, a simple scaling 
relation of the kind H,{N,T) = H,{N''T^) is found in 
the d-dimensional lattices simulated as well as in regular 
graphs. The value of a = 0.95(2) holds universally, and 
it indicates that the number of spins involved is more or 
less an extensive quantity. The temperature exponent b 
seems to depend slightly on the topology; it has a slightly 
different value for very low temperatures such as T = 0.3. 
In short, near-perfect data collapse does not seem possi- 
ble for the whole range of low temperatures simulated. 

We must emphasize that data for large values of the 
valley index must be interpreted cautiously when scaling 
the quantities in Figs. Since some runs only reach 
a few valleys, averages for larger valley indices become 
biased. This seems particularly important when looking 
at Hamming distances. Thus, we only plot data for val- 
leys reached by at least 95% of all runs. By lowering this 
threshold to, say, 80%, one will clearly see that the extra 
data points lie below the master curve. 

The strong temperature dependence of H in Fig. |31 is 
remarkable considering that the energies of the states in- 
volved are virtually independent of temperature, as im- 
plied by Fig. 2] and as one would expect for actual min- 
ima. It follows that the observed minima are nearly de- 
generate, which is also expected in spin glasses. The 
T dependence of H is likely due to the fact that at 
higher T the barriers overcome are higher, and hence 
the number m of spins involved in a quake is larger. In 
the simplest scenario where the 'downhill' part of the 
quake does not substantially change m, one can assume 
m (X H (X 7V0.95_ rpj^g energy of the barrier state, which 
scales as i? oc iV°-^*, must be carried by these m spins. 
With B cx oc iVO-^^^, we find x « 0.38/0.95 = 0.4, 
in agreement with the considerations in the previous sec- 
tion. 

The energy difference between the state of lowest 
energy in the first and the i'th valley is plotted in Fig.0] 
for 2d and 3d lattices. Ai{N,T) = Ai^N'^T'') produces 
a nice data collapse, with a = 0.87(3) and b = 0.2(1). 
Again, these exponents are found in all topologies simu- 
lated, as emphasized by the 2d data in the insert. The 
2d data bends slightly more than the 3d data in the main 
panel, preventing the ultimate data collapse. This is true 
in general of Figs. although the same scaling form 
holds across 2d, 3d, Ad, and fc-regular random graphs, 
they have slightly different curvatures. 

The linear trend in Fig. ^ means that the lowest en- 
ergy decreases (almost) logarithmically with the age, a 
feature already implicit in the early investigations by 
Grest et al. [4y|. It also implies that the energy differ- 
ence between neighboring valleys remains approximately 
constant along a trajectory. This means that the energy 
'gain' induced by a quake per participating spin is of or- 
der TV^-^^/TV"-^^ « 0(1) for the range of system sizes 
investigated. In other words, the height of the barrier 
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FIG. 4: The average difference A between the states of lowest 
energy found in the first and the i'th valley is plotted versus i. 
The data for 3d systems is presented in the main panel, where 
the combinations of L and T shown are as in Fig. 01 The insert 
shows that the same scaling exponents work equally well in 
2d, where combinations of T = 0.4 and L = 30, 100 as well 
as L = 64 and T = 0.2, 0.7 are plotted. 



which must be overcome to enter a new valley as well as 
the overlap between them depend strongly on the tem- 
perature, while the amount of energy gained is nearly 
independent of T. Hence, lowering the temperature only 
slows down the process of finding similar valleys. Fi- 
nally, we mention that the starting point Ei when mea- 
suring Ai is only slightly T-dependent. If this was not 
the case, we could not claim the similarity of valleys as 
stated above. 

To sum up, the plots of Figs. EJSl t^U us that bar- 
riers, Hamming distances, and energies of the minima 
of the valleys can all be scaled with respect to size and 
temperature in a simple way. Furthermore, the scaling 
exponents are universal; the same values of the latter are 
found in d-dimensional euclidean lattices as well as in 
regular random graphs. 



C. Residence time distribution and superaging 

The distribution of time spent in 'traps' or valleys of 
the energy landscape has, to the best of our knowledge, 
never been measured by others in simulations of spin 
glass models. The assumed form of this residence time 
distribution enters heuristic scenarios of spin glass relax- 
ation '4?, '48^1, as well as the log-Poisson description of 
non-equilibrium relaxation 1^ .38] . 

Thus, we consider the probability R{t \ t^,) that the 
residence time tr in a valley entered at age i^, be less 
than t. As with the scaling plots in the previous sec- 
tions, a clear picture emerges which lends support to the 
usefulness of the method for identifying valleys described 
in Section ITU 

Before we discuss the numerical results in any de- 
tail, we present some theoretical remarks on its ex- 
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FIG. 5: The main panel shows the probabihty that the res- 
idence time tr is less than t for valleys entered at for a 
broad range of t^,. The abscissa is the scaled variable t/f^. 
The insert shows that the pure aging t/tw scaling predicted by 
log-Poisson statistics is a fair but far less satisfactory approx- 
imation. The system is a spin glass on a fc-regular random 
graph with fc = 6. 



pected form for a log-Poisson process. Assuming a pure 
log-Poisson process for the 'log waiting time' 

\n{tr + iiu) — ln(t„) is exponentially distributed with av- 
erage 1/a, where a — dnv/d\nt is the logarithmic slope 
of the curves in Fig. ^ For any x > 0, the probability of 
tr > tw{e^ — 1) is therefore e~"^. Taking t = — 1, the 
age-scaled residence time probability distribution is 



Rit I U 



(7) 



This results also follows directly from Eq. (O by no ting 
that R{t \ tio) = 1 - Po{tw,tni + t). Since a > 1 [ij 
the distribution has the finite average {tr) — tw/{a — 
1). Equation {Tj) has been shown to fit the empirical 
residence remar kably well O. 

The t/t^ scaling form is often called pure aging, as op- 
posed to super- or sub-aging where the scaled variable is 
t/tw^ with ^ > 1 or /i < 1, respectively. Beside Fig. El su- 
peraging of the thermal correlation C{tw, t^+t) has been 
observed in numerical investigations of spin glasses ^| , 
while most experimental data show sub-aging jSCt l5lj . 
Very recent experimental work t37| shows that subaging 
is a consequence of the finiteness of the cooling rate. In 
the limit of a large cooling rate /i approaches one and 
pure aging is obtained. 

The empirical data of Fig. [S] are for random graphs of 
N = 16000 spins at T = 0.5. As with all the scaling plots 
presented so far, the quality of the fit is equally convinc- 
ing in d-dimensional lattices. To study the age depen- 
dence of R{t I tu;), the valleys entered are analyzed sepa- 
rately for a broad range of waiting time t^ = 10^, 10^. 
As the insert shows, a t/tw scaling is a fair approxima- 
tion, neglecting the small but systematic drift to the right 
as tw grows. 
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FIG. 6: The correlation between the valleys m and m + k for 
m = 5 in 4d euclidean lattices of size A'^ = 8*. The decay is 
clearly exponential for low T. Again, this is observed in all 
topologies tested. 



The main panel reveals that an excellent data collapse 
can be obtained with 



Rit I t„) = Rit/O, 



(8) 



where fi = 1.055. The growth of the residence time with 
shows that the valleys explored become more stable as 
the system grows older, which concurs with the growth 
law for the barriers given in Fig.|21 Physically, this means 
that, compared to the idealized log-Poisson case, the res- 
idence time in a valley grows even faster than the age. 
This deviation from log-time homogeneity is likely due 
to the already mentioned fact that it becomes relatively 
harder to find a new valley as the ground state is ap- 
proached [33 |. 

Allowing n to depend on T, the scaling works well 
for all temperatures up to T = 0.8. dji/dT > is ob- 
served |3^, i.e. the superaging effect becomes more pro- 
nounced as T increases. Noting that simulations slightly 
below the critical temperature, i.e. at the high end of 
the low temperature phase, are able to explore lower en- 
ergy regions in the same span of time, this strengthens 
our hypothesis that deviations from t/tw scaling are a re- 
verberation of the finiteness of the ground state energy. 
As such, these deviations can be expected to become less 
important the larger the system is, as they indeed do [ssf . 



D. Correlations 

We have investigated the overlap between configura- 
tions m and m + k for m < 5 and fc < 12 and found that 
the form 

c(m,m + fc) = CooM + (1 - Coo(m))e-'y(")'= (9) 

fits the data for our range of low temperatures in eu- 
clidean lattices as well as in random regular graphs. Only 
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in 2d do we find significant deviations, as expected con- 
sidering that the aging is interrupted in this case ^3 • 
an example, Fig.Elshows c(5, 5 + fc) in Ad models. A simi- 
lar plot for c(l, l + fc) in Sd systems can be found in Fig. 4 
of Ref. ^3 ) where the connection to the experimentally 
available non-equilibrium exponent A is verified. 

If the limiting value for fc — > cxd of c{m,m + k), 
Coo(w), were independent of m, it would coincide with 
the Ed wards- Anderson order parameter. Empirically, we 
find a small m dependence, for which we have no physi- 
cal interpretation. The exponent 7 also has a small and 
almost linear m dependence. As Eq. @ remains a heuris- 
tic approximation of limited value, we have not pursued 
the size and temperature dependence of the parameters 
involved in Eq. ©. StiU, we find it interesting that a sim- 
ple exponential parametrization in the number of quakes 
k accurately describes the data. 



E. Hitting time for the minima 

Having argued that our method of partitioning the 
sampled states into valleys leads to consistent results and 
provides useful insight into the coarse structure of com- 
plex energy landscapes, it is natural to take a first look at 
the rich internal structure of the valleys. The existence 
of such a structure is implied by the wide distribution of 
residence times, which requires a matching distribution 
of internal energy barriers. Direct evidence is presented 
in this section, where we consider the 'hitting time' thu 
elapsing from the time of entry to the time where the 
state of lowest is encountered. 

Consider r = thu/tr, the fraction of the residence time 
spent 'searching' for the global minimum. We expect r to 
be close to zero in a structureless valley, where the global 
minimum is reached soon after entry, and close to one in 
the opposite limit of a rugged valley with many internal 
minima. In measuring the distribution of r, we keep track 
of the age t^ at the time of entry, by selecting entry times 
for a range of different tw From the outset, the empirical 
distribution G{t \ t^) could depend on i^, but, as shown 
by Fig. there is practically no dependence, except for 
T close to the critical temperature. 

The right-skewed form of G for lower T values (the 
picture is the same for any T < 0.6) implies that by far 
the greatest part of tr is spent before hitting the low- 
est minimum: the system explores several sub-valleys, 
each identified by its own local minimum, while it slides 
towards states of lower energy. Equilibrium-like fluctua- 
tions can only occur after the state of lowest energy which 
remains at the bottom of the valley is hit. This time in- 
terval is small compared to the total residence time, but 
nevertheless grows steadily with the valley index. 

Since the measured shape of G at low T is almost 
invariant, the distributions of both tr and thu must both 
scale as tf^ . We conclude that the process of hitting the 
lowest minimum in a valley at low T is limited by inter- 
nal barrier crossing events, and that these barriers them- 
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FIG. 7: The location of the energy minimum within a valley 
entered at time t„ in 3d spin glasses with — 16^ spins. The 
abscissa is the time relative to the entry and exit times. The 
figure illustrates the highly right-skewed distribution found 
for all low T, independent of t^. The insert shows that for 
higher T the distribution flattens out as the age tw of the 
system grows. 

selves grow as the age increases and the valleys become 
larger. The landscape geometry is thus invariant under a 
time dilation, and hence self-similar, in agreement with 
the general properties of a log-Poisson process. The tem- 
perature independence of the distribution at low T leads 
to the same conclusion: Lowering the temperature means 
exploring smaller valleys, which nevertheless retain the 
same internal structure within a range of low tempera- 
tures. 

As the temperature gets closer to the critical temper- 
ature, a waiting time dependence appears: G(r) remains 
peaked around r ~ 0.9, but smaller values of r become 
gradually more probable. The corresponding flattening 
of the distribution of G(t) for T — 0.8, a high tempera- 
ture, is shown in the insert of Fig. |7| Since the number 
of potentially new valleys diminishes in time because of 
the lack of new low energy records, G must eventually 
become highly left-skewed. While only a very slow trend 
towards this situation can be seen in the main panel, the 
data for the higher temperature in the insert seems to be 
heading in that direction, in accordance with equilibra- 
tion happening much faster at high T than low T. 

We end this section by noting that, like most other 
quantities presented in this paper, similar results for 
G{t I tw) are found across all systems simulated. 



VI. COARSE GRAINED LANDSCAPE 
DESCRIPTION 

Having argued that log-Poisson statistics is a (slightly) 
idealized description which leads to pure aging and other 
dynamical features found in experiments, we turn to the 
physical mechanism behind the selection of the attrac- 
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tors and, more broadly, consider the implications of our 
results for pertinent mesoscopic models of complex land- 
scapes. We do not attempt to deal with domain growth 
and other real space issues in any detail. These aspects, 
though very important for a complete understanding of 
complex relaxation, are too strongly connected to quasi- 
equilibrium issues which lie beyond the scope of this pa- 
per. 

The applicability of log-Poisson statistics to barrier 
records only implies that the trajectory fully decorrelates 
between successive events, as can be expected for an ac- 
tivated process in a landscape with many local minima. 
The new and important information about the landscape 
geometry lies in the fact that barrier records are required 
in order to find new minima records. This explains why 
the log-Poisson statistics is equally relevant for barriers 
as for minima. Secondly, and more importantly, it im- 
plies that the least stable, or marginally stable among the 
available attractors are those selected. Indeed, if at some 
stage, very deep minima were to follow a shallow bar- 
rier, subsequent barrier records would not produce new 
states of record-low energy, and the log-Poisson descrip- 
tion would fail. 

The close match between the depth of a valley and the 
magnitude of the barrier record giving access to it im- 
plies a quasi-continuum of available attractors and was 
first noticed in connection with the noisy dynamics of 
a driven dissipative system [T^ . where the phenomenon 
was dubbed noise adaptation. That complex memory be- 
havior is linked to marginal stability of metastable attrac- 
tors has long been known for noiseless models of Charge 
Density Waves The extremely simple automaton 

model of Tang et al. 52] describes a sheet of elastically 
coupled 'balls' driven along a sinusoidal potential by a 
pulsed external force. A recent study of a noisy version 
of this model 10] shows that its non-equilibrium dynam- 
ics is described by log-Poisson statistics and that the age 
of the system can be reset by a change of the elastic cou- 
pling constant. 

The physical origin of the bias towards shallow attrac- 
tors in thermal glassy dynamics is likely to be entropic, 
i.e. simply the fact that shallow attractors vastly out- 
number deeper ones, in line with the general observation 
that quenches usually produce poor minima. To sup- 
port such bias for a range of low temperatures, the den- 
sity of energy minima must dwarf the Boltzmann factor 
and hence increase at least exponentially with the energy. 
This concurs with the outcome of numerical exhaustive 
investigations of the local configuration space structure 
of different glassy systems ^23,, 24, 41, 43, 44, 53, 54] 
which were performed with the lid method |23L |42| . In 
all cases the local density of states and the local density 
of minima are nearly exponential functions of the energy. 

An exponential density of states in connection with 
activated dynamics implies a dynamical glass transition. 
This exact feature was built into the tree model of com- 
plex relaxation proposed and analyzed by Grossmann et 
al. [E^ and later studied in more detail in |5^ [E^]. It 



is also incorporated in the even simpler trap model of 
Bouchaud [43 which, nonetheless, is very different from 
tree models in one important respect: trees have the low- 
est connectivity possible for a connected set, while each 
trap of the trap model is connected to all others. 

Log-Poisson statistics only applies as long as new and 
gradually more stable valleys remain available to the dy- 
namics. The simplest way of modeling inequivalent val- 
leys is through a hierarchy of energy barriers separat- 
ing degenerate states, which are organized in either a 
linear array or in a tree graph: Beside the mod- 
els alr eady mentioned , th e latter approach is followed 
in m Many important features of 

complex relaxation can be reproduced in this model, but 
not the fact that in many systems, including spin glasses, 
the energy decreases logarithmically with the age. This 
can be achieved by introducing non-degenerate local min- 
ima, as done in the so-called LS tree 7, 63]. The minima 
have energies which on average decrease linearly with the 
size of the barrier overcome, i.e. in the same overall fash- 
ion as Fig. 01 

In Bouchauds model, trap energies are exponentially 
distributed and hence non-degenerate. That deeper min- 
ima are gradually explored is a statistical consequence 
of the (assumed) infinite average of the residence time 
in a trap. If traps and valleys can be identified, the 
results of Section IV CI are at variance with this inter- 
pretation. The average residence time, which equals tw 
for a = 2, is, in practice, slightly lower than 1^. This 
is again reminiscent of the situation encountered in tree 
models Hillll IS El HI El . The distribution of bar- 
riers in tree models has a lower cut-off, unlike the fractal 
description of configuration space of Dotsenko 64, 65|, 
which better captures the growing importance of gradu- 
ally smaller barriers as the temperature decreases. 

A last important issue is the connection between the 
energy barrier separating two configurations and the dis- 
tance between them. For spin glasses the relevant met- 
ric is the Hamming distance, which, according to Figs.|2| 
and^ on average bears a linear relationship to the energy 
barrier. A similar result was found in numerical work 
on the SK model by [Hil], and by the lid-method 
(i.e. exact exhaustive enumeration) in jlj for short range 
spin glasses. Since the latter investigations deal with the 
small scale structures inside a 'pocket', while the present 
ones are concerned with the large scale structures ex- 
plored by the non-equilibrium dynamics, the agreement 
in their outcome is further evidence of a self-similar land- 
scape structure. Interestingly, the largest distance which 
can be achieved for a fixed energy barrier grows exponen- 
tially with the barrier |2^ . A linear relationship between 
energy and Hamming distance is assumed in tree models 
of aging dynamics, see e.g. [^J, and also in the barrier 
model 58]. This is, however, not a crucial assumption 
for aging, and other types of functional relationships can 
also be utilized ,3] successfully. 
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VII. SUMMARY AND OUTLOOK 

In this paper we have presented a general 'non-invasive' 
statistical method for complex energy landscape explo- 
ration, especially designed to provide information on the 
metastable states actually explored by an unperturbed 
aging process following a quench. The method has been 
thoroughly tested on Ising spin glasses, and the results 
obtained both match and extend the established knowl- 
edge about spin glasses. In particular, most quantities 
investigated obey simple scaling laws with universal scal- 
ing exponents. The view which emerges is that non- 
equilibrium aging dynamics is steered by energy barrier 
records, which are the only events capable of opening the 
route to new valleys. Since this description has previ- 
ously been shown to apply to driven dissipative models 
and to evolution modeling, a possible unified theory for 
non-equilibrium glassy dynamics seems within reach. 

A more complete picture can be obtained by looking 
more closely into the fluctuation dynamics which may 
differ across different systems. For spin glasses we have 
argued that real space domains of (pseudo) thermalized 
spins relax independently as long as the system as a whole 



remains in the same valley. The present method opens 
the possibility of identifying the quasi-equilibrium clus- 
ters as defined by the dynamics itself: These clusters are 
separated by a backbone of spins whose orientation re- 
mains fixed within each valley and changes slowly from 
one valley to the next, as seen in Fig. El After the low- 
est energy state in the valley has been hit and before 
the next valley is entered, there is no drift towards the 
global minimum. Hence, the spins fall into two categories 
only: those which are frozen, and those which fluctuate in 
a quasi-equilibrium fashion. The quasi-equilibrium clus- 
ters can thus be extracted and their statistical properties, 
such as e.g. the density of states, can be studied for each 
cluster separately. 
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